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ABSTRACT 

In previous work by the author, a generalization of Godunov’s method for systems of conservation laws 
has been developed and analyzed that can be applied with arbitrary time steps on arbitrary grids in one space 
dimension. Stability for arbitrary time steps is achieved by allowing waves to propagate through more than one 
mesh cell in a time step. In this paper the method is extended to second order accuracy and to a finite volume 
method in two space dimensions. This latter method is based on solving one dimensional normal and tangential 
Riemann problems at cell interfaces and again propagating waves through one or more mesh cells. By avoiding 
the usual time step restriction of explicit methods, it is possible to use reasonable time steps on irregular grids 
where the minimum cell area is much smaller than the average cell. 

Boundary conditions for the Euler equations are discussed and special attention is given to the case of a 
Cartesian grid cut by an irregular boundary. In this case small grid cells arise only near the boundary, and it is 
desirable to use a time step appropriate for the regular interior cells. Numerical results in two dimensions show 
that this can be achieved. 
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1. Introduction. 

Hyperbolic systems of conservation laws in two space dimensions, such as the Euler equations of gas 
dynamics, have the form 

*.+/0O» + «(«),= 0. (1.1) 

Here u is a vector of m variables and / and g are the flux functions. In order to solve these equations 
numerically, the domain is subdivided into computational grid cells (or intervals in one dimension) and an 
approximation to the solution over each cell is updated in each time step. With traditional explicit methods, the 
appropriate update is based on cell values in only a few neighboring cells. For example, with a three-point 
method in one dimension, only the nearest neighbor cells on either side can effect the solution. Most finite 
volume methods in two space dimensions also limit the domain of dependence to the nearest neighbors. 

This imposes an obvious limitation on the time step due to the CFL condition, which states that a 
necessary condition for convergence is that the domain of dependence of the numerical method include the 
domain of dependence of the true solution, at least in the limit as the grid is refined. For a three-point method in 
one dimension, this means that the time step k must be restricted by 

V-T^-IJmJSl (1-2) 

™miii 

where h„^ is the length of the smallest cell and s mix is the maximum wave speed arising in the problem at 
hand. The quantity v defined by (1.2) is called the Courant number. 

This time step restriction may be reasonable on a uniform grid. An analysis of the truncation error will 
often show that (1.2) is consistent with the choice of time step required for the temporal accuracy to agree with 
the spatial accuracy. On the other hand, if the grid is not uniform we would like to replace h^ in (1.2) by h ive , 
some average cell length. If « /t lve then (1.2) imposes a severe time step restriction that may greatly 
reduce the efficiency of the method. 

This situation arises, for example, in certain shock tracking procedures (e.g. [16]). In addition to a 
uniform grid, moving points are introduced that may fall arbitrarily close to the fixed points, creating small 
cells. Similar difficulties are encountered in the moving mesh method of Harten and Hyman[12], where they 
must take special care to insure the proper separation of mesh points. 

In more space dimensions, the problem of small mesh cells is difficult to avoid in several common 
situations. In particular, irregular boundaries cutting through a regular grid (as in Figures 9.1 and 9.3) result in 
cells with arbitrarily small area. We would like to choose the time step based on the regular grid cells rather 
than being artificially restricted by the boundary cells. Small cells are often avoided by using body fitted 
coordinates, so that the grid cells have more uniform areas. However, this creates additional difficulties in first 
generating and then computing with the irregular grids. For complicated multicomponent geometries, 
particularly in three dimensions, the generation of such grids may be very difficult, requiring complicated grid 
patching. This has lead several workers to reconsider the use of Cartesian grids (e.g. [6,25,28]). 

Similar difficulties arise if a shock is tracked through a fixed grid, at the boundary of a region of mesh 
refinement, or where independent grids are patched together. On general unstructured grids, such as 
triangulations of arbitrary regions, the cell areas may also vary widely. 
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In this paper a finite volume method is introduced that avoids these difficulties by dynamically extending 
the domain of dependence in such a way that the CFL condition is always satisfied. This is accomplished by 
fully exploiting the wave propagation that is inherent in hyperbolic equations. A one dimensional Riemann 
problem is solved normal to each cell interface (meaning the initial data is piecewise constant with a jump 
discontinuity at the interface). Solving Riemann problems forms the basis of Godunov’s classical method and 
numerous finite volume methods, but in most such methods the Riemann solution is used only to obtain a flux 
across the boundary. This flux is then used to update only the neighboring cell values. 

In the method proposed here, the structure of the Riemann solution is used to a much greater extent. The 
Riemann solution consists of waves (shocks, contact discontinuities, or rarefaction waves) that propagate away 
from the interface, each carrying an increment in the flow variables. As a wave propagates, it modifies the flow 
variables in the cells it passes through. By allowing it to pass through more than one cell, the CFL condition 
remains satisfied with larger time steps. This gives a generalization of Godunov’s method applicable for larger 
time steps on arbitrary grids. In one space dimension, this first order accurate method has been developed and 
analyzed in a series of papers by the author[14,15,16). This method is similar in spirit to the "transport collapse" 
method of Brenier[l,2,3]. In the present paper, the method is extended to second order accuracy and arbitrary 
grids in two space dimensions. 

In addition to propagating waves normally from cell interfaces, it is also possible to incorporate 
tangential wave propagation by solving auxiliary Riemann problems in the orthogonal direction. This leads to 
increased stability and enhanced modeling of two dimensional phenomena and propagation of shock waves that 
are not aligned with the grid. 

An additional advantage of the method proposed here is that boundary conditions can be incorporated 
into the algorithm relatively easily. Farfield nonreflecting boundary conditions and solid wall (or moving 
piston) boundary conditions for the Euler equations will be discussed here. 

The basic idea of the method is described first in Section 2 for a scalar equation in one space dimension. 
A new formulation of the Lax-Wendroff method is presented and then converted to a "high resolution" method 
(one with sharp nonoscillatory discontinuities and good nonlinear stability properties). The formulation of this 
method easily allows extension to arbitrary grids and time steps and incorporation of boundary conditions. 

In Section 3 the method is extended to systems of conservation laws in one dimension. Boundary 
conditions for the Euler equations are discussed in Section 4 and then some one-dimesional numerical results 
are presented in Section 5. 

In the remainder of the paper the method is extended to a finite volume method in two dimensions. The 
special case of a Cartesian grid with irregular boundaries receives particular attention. In Section 9, numerical 
results are presented for a shock tube problem oblique to the grid and single Mach reflection from a ramp. 

2. One dimensional scalar equations. 

We first consider the one-dimensional scalar conservation law 

«,+/(«)* = 0 . ( 2 . 1 ) 

with /(«) assumed to be a convex function. The computational grid is given by jc 0 <x 1 < • • ■ <x M . The 
mesh spacing can be arbitrary and we set A; = x 1+1 -x ( . The approximate solution on the mesh cell (x, ^t, +1 ) at 
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time t„ is denoted by t/". The time steps can also be arbitrary and varying. Since the methods are one-step 
methods and will be described over a single time step, we will denote the current time step simply by 
k = t n+1 - t n and use the short hand abbreviations (/,• = {/" and l/,- s {/" +1 . 

For each value of j , set 


Wy+ 1) ~f(Uj)) I (U j+ , - (/;) if t/y * f/y+i 
“/'(t/y) if t/y = U j+l . 


( 2 . 2 ) 


Then sy is the speed at which the discontinuity between U y and U y + i propagates according to the Rankine- 
Hugoniot relation. For convenience at this point, suppose that Sj> 0 for all j . Then we set 


Uj(X,t n ) 


U j for x <Xj 

' Uj + (x -x j+V2 )Oj for x j <x < Xj+i 
Uj+l for * Sxy +1 


(2.3a) 


where Xy +1/2 = y(xy +x J+1 ) and Cj is the slope of Uy over the j th interval, as yet unspecified (see Figure 2.1). 

Based on these functions, we define an approximate solution u(x ,t ) over the time interval < t < t n+1 . At time 
t n , u(x ,t n ) is simply the piecewise linear function that agrees with Uj(x ,t n ) on the j th cell: 


U (X ,t n ) = Uj(x,t n ) for Xj<X <Xy +1 


(2.4a) 


For t >t H we let the wave form (2.3a) propagate with speed s jy defining 


Uj(X,t)=Uj(X -Sj(t-t H ),t n ). 

These waves are then combined via linear superposition to obtain the approximate solution 


(2.3b) 


u(x,t) = u(x,t„) + 2 ( Uj(x,t ) - Uj(x,t H )). (2.4b) 

j 

The sum in (2.4b) is over all j but at any given point only a finite number of terms will be nonzero since 
Uj(x,t H ) is constant outside of (Xjjc j+l ). 

The approximation U, at time /„ +] is obtained by projecting the approximate solution u(x onto the 
grid, as in Godunov’s method (i.e. by averaging over the mesh cell (x,-,x i+1 )): 


— 1 

= u(x,t n+1 )dx 

= U i + -r'LC (Uj(xj H+l )-Uj(x,t H ))dx. (2.5) 

«/ y * 

Note that t/, is obtained from f/j by adding to it contributions from several neighboring cells. Figure 2.1 shows 
the typical contribution to f/,-, indicated by the shaded area. 

For a linear problem u, + m, = 0 with constant speed s, the approximation u(x,t) is in fact the true 
solution with initial data u(x,t„). For nonlinear problems the waveform (2.3a) should be distorted for t > f„, 
and so (2.3b) represents a local linearization of the flux function. The superposition (2.4b) is a further 
linearization of the nonlinear problem. In spite of this, it is possible to achieve second order accuracy with 
appropriate choices of Oj. 
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There are various ways to specify the slope Oj. Three possibilities are: 


II 

o 

(Godunov) 

(2.6a) 

Vj W = (Vj+ 1 - Uj) / (Xj+3,2 - X >+1/2 ) 

(Lax-Wendroff) 

(2.6b) 

of = minmod (cj 4 w , off) 

(minmod) 

(2.6c) 

where the minmod function is defined by 




minmod (a Jb ) = ~ (sgn(a ) + sgn(b )) min( I a I , I b I ). 


The names assigned to the first two choices of o are based on the fact that the method reduces to one of 
these standard methods on a uniform grid with Courant number v < 1. Then each [/, is influenced by at most 
two waves and the method can be rewritten as a three-point scheme. Setting v ( - = s,jfc/h and computing the 
integrals in (2.5) explicitly, we obtain (still assuming Sj > 0 for all j) 

Ui = £/, - v,_i(f/i - + ~V;_i(l - v, _!)©,_!* - -jv,(l - v.taft. (2.7) 

This can be rewritten in conservation form as 

Ui = U i -j ; (F i -F i _ l ) (2.8) 

where 

Fj =f(Uj) + y(l - Vj)Sj(Sjh. (2.9) 

If Cj = 0 this gives the first order upwind scheme, which agrees with Godunov’s method for the problem 
considered here. With slope (2.6b), this is easily seen to be the standard Lax-Wendroff method. 

The Lax-Wendroff method is second order accurate on smooth solutions, but it generally develops 
oscillations and possibly instabilities near discontinuities or sharp gradients. The minmod choice of slopes 
(2.6c) is designed to avoid this difficulty. Van Leer[27] used a similar geometric approach with limited slopes 
to define his second order accurate MUSCL schemes. Later Goodman and LeVeque[10,18] used the minmod 
slope (2.6c) to derive a total variation diminishing (TVD) version of van Leer’s method. The method (2.7) with 
this choice of a is closely related to these methods and to a variety of other "high resolution" schemes that have 
recently been proposed (e.g. [1 1,19,20,29]). In particular, for a linear problem with f(u) = au and a > 0, (2.7) 
with minmod slopes is precisely the minmod flux-limiter method introduced by Roe (see Sweby[26]) and hence 
is TVD. 

It is interesting to note that the minmod method proposed here is not TVD on nonlinear problems. It is 
possible to produce examples where the variation increases slightly, although detrimental oscillations or 
instabilities have never been observed. Indeed, the lack of TVD may be an advantage since any TVD method 
must reduce to first order accuracy at extreme points[19]. In numerical experiments, the present method 
appears to perform better at such points. 

The main advantage of this approach, however, is that we can easily go beyond formula (2.7) and apply 
the method as described by (2.5) with arbitrary grids and time steps. 
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Before discussing systems of conservation laws, we point out a slightly different formulation of the 
method that helps simplify its implementation and extension to two space dimensions. We can decouple the 
wave propagation into a two stage procedure. In the first stage we propagate the wave u} l \x,t) given by (2.3) 
with Oj = 0, corresponding to the first order algorithm. In the second stage we propagate a second order 
correction wave. 


U/ 2) (X,I) = Uj{X,t) - u/'W) 

with the same speed Sj. This wave has integral zero and initially has the form 


uf 2 \x,t n ) = 


(x -Xj +ia )Gj 
0 


for xj <x <Xj + i 
otherwise 


When viewed in this way it is clear how the second order correction relates to the flux corrected transport 
(FCT) algorithm of Boris and Book[3] and Zalesak[30] and to the ideas of Roe[22]. The results of the first 
order accurate scheme are corrected by redistributing some of the flux between grid cells. 


In the above discussion we assumed that Sj > 0. If Sj < 0 we must replace (2.3) by 


(x i 


Uj 


for x < ; tj+i 


Uj + (x - x J+3a )<Jj for x j+ i <x < Xj +2 


( 2 . 10 ) 


U 


/+! 


for jc >*y + 2 


in order to obtain a method that reduces to the Lax-Wendroff or minmod method in the special cases described 
above. In the latter case we also redefine 


of = minmod (of™, of™). 

Note that the linear segment in Uj (jc is always behind the discontinuity as it propagates. 


3. Systems of conservation laws in one dimension* 

We now develop the method more generally for hyperbolic systems of conservation laws (2.1) with 
u € R m . We assume that the Jacobian matrix f\u) has real eigenvalues ^(u) for all u and that each 
characteristic field is either genuinely nonlinear or linearly degenerate. Then the Riemann problem, consisting 
of (2.1) with piecewise constant initial conditions 


(Uj x<x J+l 

«(*,'„)= | ^ +i x > Xj+i (3-D 

has a solution for //f/ /+] - {^//’sufficiently small. This similarity solution consists of m + 1 constant states 
separated by shocks, contact discontinuities, or rarefaction waves. (See Lax[13] for a description of the 
theory.) 

Although the method can be adapted to use this exact Riemann solution, we will base the method on an 
approximate Riemann solution of the type introduced by Roe[21]. The Riemann solution of (2.1) with data 
(3.1) is replaced by the solution of the linear problem 


u, +AjU x = 0 


(3.2) 
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where A; =Aj(Uj,U J+1 ) is an appropriate average value of f'(u). This gives a good approximation, 
particularly for smooth solutions where flU y +l - Uj// = 0(h). Denote the solution of this problem by uf l \x ,t ) 
for t >t n . Then u.p\x,t) consists of m + 1 constant states separated by discontinuities travelling with speeds 


which are eigenvalues of Aj. The jump in 


s n <Sj 2 < ■■■ <s jm (3.3) 

uj l) across the p th such wave is an eigenvector Rj p of A so that 


AjRjp — s jpRjp (3.4) 

and V y + i - Uj is decomposed as 

U J« ~ u j = l. V (3-5) 

P =i 


To obtain high resolution we will also introduce a correction wave uj 2 \x,t) by generalizing the scalar 
case discussed in Section 2. We break «/ 2) up into m correction waves uj- p \p = 1,2, ••• , m each of which 
travels with speed s jp and has an initial form that is piecewise linear and nonzero over a single mesh interval. 
The choice of interval again depends on the sign of s jp : 


For Sj p > 0: 


ufp\x,t n ) = 


f(x -x j+V2 )o Jp 

10 


for Xj< x <Xj +1 

otherwise 


(3.6a) 


f(j 

For s j P ^ 0: u J £ ) (x,t n )= K 


(* - Xj+ s/z)Gj P for x j+l <x< x j+2 
otherwise 


where a jp is now a slope vector with m components, to be defined below. We now set 

u J a \x,t)= £ ug\x,t)= £ ugXx -s Jp (t -t n ),t H ) 
p = i p=i 

and 


(3.6b) 


(3.7) 


Uj(x,t) = Uj m (xj) + uf 2 ) (x,t). 


(3.8) 


Based on these functions Uj(x,t), we can define an approximate solution u(x,t) just as in the scalar case via 
(2.4). Averaging u(x ,t n+l ) over the i th interval as in (2.5) gives the numerical solution U, at time t n+l . 

If s 0 then «y®(r))=0 for all j. In this case the method reduces to the large time step 
generalization of Godunov’s method (with Roe’s approximate Riemann solver) that has been previously 
studied([14, 15,16]). This method is conservative for any mesh ratio as is shown in [14] where the method is 
rewritten in conservation form. Note that the method remains conservative with arbitrary slope vectors since 

j“ u} 2 \x,t) dx = 0 for any a jp . We can view the function Uj(x,t) for t > t n as a new approximate Riemann 

— oo 

solver which gives second order accuracy for an appropriate choice of a Jp . Because of this the high resolution 
method with nonzero a ip can also be written in conservation form following Section 4 of [14], 


Corresponding to the three slope choices (2.6), we now have 




(Godunov-Roe), 


(3.9a) 
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ajp W = */p 1 (*/+ 3/2 - ^>+ 1 / 2 ) (Lax-Wendroff), (3.9b) 

a" = minmod (a£ w , a}*) (minmod), (3.9c) 

where in (3.9c) J =j - sgn(sj p ) and the minmod function is applied to vectors componentwise. Either of the 
last two choices appears to give a second order accuracy on smooth solutions. In order to achieve stability and 
nonoscillatory behavior near discontinuities, the minmod slope (3.9c) is recommended and has been used in the 
numerical experiments reported below. 

This method appears difficult to implement from the description just given. In fact, it is quite simple if 
we consider each wave arising from each Riemann problem in turn and compute its effect on the solution. Due 
to the linearization of wave interactions that is used here, each wave can be propagated and averaged onto the 
grid independently of all other waves. For the case a jp =0 this algorithm has already been described in some 
detail in [15] and [16]. This gives propagation of the waves uj p \ 

The correction waves ujj® simply give a redistribution between mesh cells. In particular, if the mesh is 
uniform, then the nonzero portion of uj^ overlaps at most two cells, independent of the Courant number, and 
the second order correction is easily implemented. 

For linear systems of equations, solution of the Riemann problems gives resolution into the characteristic 
variables and u(x,t) is again the exact solution with the piecewise linear data u(x,t n ), for all t>t n . 
Consequently, the method gives second order accuracy with any size time step in this case. For nonlinear 
problems with smooth solutions, the superposition (2.4) can be shown to be consistent with second order 
accuracy. Towards this end, a semidiscrete version of this superposition (discretized only in time) is 
investigated in [17] and shown to produce an approximate solution with 0 (£ 3 ) error over a time step of length 
k. 

4. Boundary conditions for the Euler equations. 

One advantage of the approach used here is the ease with which boundary conditions can be derived and 
implemented, based on the physics of the problem. This will be illustrated for the Euler equations at a free field 
(non-reflecting) boundary and at a solid wall (or moving piston) boundary. 

The Euler equations of gas dynamics have the form 

p e v 

pv + p v 2 + p =0 (4.1) 

E Jr P +P ^ V \x 

where p, v,p, and E represent the density, velocity, pressure and energy, respectively. These quantities are 
related by E = pl(y- 1) + ypv 2 with y= 1.4 for air. 

Non-reflecting boundary conditions are easily achieved by simply allowing all waves which cross the 
boundary to pass out of the region, with no incoming waves. In a subsonic problem these boundary conditions 
are not stricdy correct. For example, if two shocks propagate out of the domain at different times, they may 
interact (in the exterior) and give rise to a weak rarefaction wave which should reenter the domain at a later 
time. Correctly modeling such behavior requires boundary conditions which are nonlocal in time. Simply 
letting the waves pass out appears to be the best that can be achieved with local boundary conditions and is 
easily implemented with arbitrary grids and time steps. 
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At a solid wall boundary the correct boundary condition is v = 0. This gives rise to a boundary Riemann 
problem which can be solved for the incoming wave. Note that the boundary is characteristic in this case, so 
that there is only one incoming wave. 

Another way to impose this boundary condition for the PDE is to extend the domain beyond the boundary 
and specify initial conditions obtained by suitably reflecting the values from the interior. For example, if we 
wish to compute over time k and s mtx represents the maximum wave speed, then it suffices to extend the 
domain from x = 0, say, to x* = -ks mtx and set 

pC*) = p(-*) 

v(x) = -v(-x) (4.2) 

P(x) = p(-x) 

for x m <x < 0. Solving the problem over this extended domain (using, for example, non-reflecting boundary 
conditions at x = x * ), will give the desired solution for x > 0 with the boundary condition v (0,/) = 0 satisfied. 

This latter approach is easy to implement numerically on arbitrary grids. New grid points are introduced 
exterior to the original domain by reflecting the interior points which fall within distance ks mn of the boundary. 
The values on these mesh cells are obtained by reflecting the interior values as in (4.2). Then any wave leaving 
the original domain during the time step is matched by a "reflected" wave which enters at the same point in 
time. The jump in v in the entering wave is the negative of the jump in v in the exiting wave (due to the 
symmetry in solutions to the Riemann problem under the reflection (4.2)) and hence the boundary condition is 
always satisfied. This same symmetry exists in Roe’s approximate Riemann solver and the slopes o JP will also 
be symmetric, so that the boundary condition remains satisfied for the method discussed here. Again the 
procedure is easily implemented for arbitrary time steps and meshes. Finally, it is clear that the method has the 
same order of accuracy near the boundary as in the interior, since after the reflection the boundary essentially 
disappears and we are simply using the interior scheme everywhere. 

If the solid wall boundary is moving (as for example in a piston problem), then boundary conditions are 
also easily derived. If the piston is located at x = z„ at time f„ and is moving with speed s 0 (assumed constant) 
for t„ < t < t B+1 , then the mesh cells for z„ < x < z„ + s mn k are reflected to the left of z n with grid values 

p(z B -x) = p(z„+x) 

v(z»-x) = 2so-v(z„ + x) (4.3) 

P(.z*-x)=P(z n +x) 

The waves generated from these cell boundaries will interact with the interior waves in such a way that 
v (z„ + s 0 (f - f„),f ) = 0 will be satisfied for t H < t < t H+l . 

5. Numerical results in one dimension. 

On a uniform grid with v < 1 the present minmod method performs much like the minmod flux limiter 
method discussed by Sweby[26]. Figure 5.1 shows some results for a shock tube problem with the Euler 
equations and initial data with v = 0 everywhere and p =p = 1 for x < 0.4, p = p = 3 for x > 0.4. The solution 
at three different times is shown. Note that the shock reflects off a solid wall boundary at x = 0 and then 
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interacts with the contact discontinuity. Nonreflecting boundary conditions were used at x = 1, where the 
rarefaction wave begins to leave the region. 

Figure 5.1a shows the density p computed using the first order accurate method with Godunov slopes 
Cj p = 0. Figure 5.1b shows the results with minmod slopes. In each case h = 0.02 and k = A/2 corresponding 
to v = .75. The solid line in each figure is the solution computed with 201 mesh points using minmod slopes. 

Figure 5.2 demonstrates the ability of this method to deal with nonuniform meshes. In this problem 
smooth initial conditions were specified: 

p(jc ,0) = 1. + 0.1cos(2.5jc), v ( x ,0) = 0.1sin(5x), p (. x ,0) = 1.8 + 0.1cos(2.5x), 

again with a solid wall boundary at x = 0 and nonreflecting boundary conditions at x = 1. The solid line shows 
the resulting density at t = 0.1 computed with 51 equally spaced mesh points and k - h = 0.02, giving v = 1. 
The squares show the solution computed with 51 mesh points Xj = (0 — l)/i) 3 densely clustered near x = 0. 
Again k =0.02 but now the Courant number, based on the smallest mesh width h u is roughly 2500. The 
numerical solution does not appear to be adversely affected by this. 

6. Finite volume methods in two dimensions. 

In two space dimensions, the wave propagation approach can be used to rewrite a standard finite volume 
method in a way that can be easily extended to allow larger timesteps than are normally possible. As noted in 
the introduction, this might be valuable with nonuniform grids where some of the grid cells are considerably 
smaller than average. Standard explicit methods have a stability restriction that limits the time step in terms of 
the smallest cell area. Of particular interest is a uniform Cartesian grid in which some of the grid cells are cut 
off by an irregular boundary, a tracked shock, or an area of mesh refinement. To some extend the discussion 
will be specialized to this situation, but the ideas can be applied on more general unstructured grids as well. 

In one space dimension, we have seen that nonuniform grid cells are easily handled by allowing waves to 
propagate entirely through a small cell and out the other side. The same approach can be taken in two 
dimensions with waves arising at cell interfaces. The basic algorithm will be described on a general 
unstructured grid. 

Associated with each grid cell C } is an approximation Uj to the solution u(x,yj n ) over this grid cell. 
We also define Aj to be the area of Cj and hy to be the length of the face separating adjacent cells C x and C r 
We begin by describing a Godunov-type scheme and will later introduce high resolution corrections based on 
piecewise linear approximations. 

In a finite volume method, fluxes across the cell faces are defined and used to update the cell values on 
either side of the face. Consider, for example, the face separating C, and Cj in Figure 6.1. A flux can be 
defined by solving a one-dimensional Riemann problem normal to the face. Define new variables £ (normal to 
the face) and tj (tangential) by 

£ = ax + Py, rj = -px+ay (6.1) 

with a = cos0, p = sin0 where 0 is the angle of the cell interface. The conservation laws 


u,+f(u) x + g(u) y = 0 


( 6 . 2 ) 
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are transformed to 


u, + F(u)f. + G(u) 1) = 0 

where 


(6.3) 


F(u) = af(u) + Pg(u) (6.4) 

G(u) = ~Pf(u) + ag(u). 

In the new coordinates u is constant in q and so (6.3) reduces to a one-dimensional Riemann problem 

u,+F(u\ = 0 (6.5) 

with left and right states Uj and Uj. For rotationally invariant equations such as the Euler equations this is 
particularly simple since, after a change of dependent variables to rotate the velocity field, the form of F(u) 
agrees with / (u) and hence a single Riemann solver suffices for all angles 0. 

As in the one-dimensional algorithm, we use Roe’s approximate Riemann solver and obtain a set of 
waves traveling with speeds s lt s 2 , • • • ,s n in the ^-direction. As before we denote the jumps in u across these 
waves by the vectors R j, • • • ,/? m , so that 


U J -U i ='ZR p . 

p 


( 6 . 6 ) 


Figure 6.1 shows these waves propagating from the interface in a situation where the waves lie entirely within 
neighboring cells. 


We think of u being incremented by -sgn (s p )R p at each point swept out by the p\h wave. Consequently, 
the average value t/ ( of u over the cell C, in Figure 6.1, for example, will be incremented by 


\s l \kh ij 

A, 


R i- 


(6.7) 


The factor multiplying R t is the ratio of the area swept out by the 1-wave to the total cell area. Similarly, Uj is 
incremented by 


sjchj. 


W + 


sjchjj 

Ai 


H? 3 ). 


( 6 . 8 ) 


These increments to (7, and Uj can be interpreted more traditionally in terms of a flux-difference 
splitting. The total increment to u (which has been split between C, and C y and averaged over these cells) is 
the sum of A; times (6.7) plus A j times (6.8), giving 

-khijUfti + S2^2 + ^3^3]- (6.9) 

A basic property that Roe’s approximate Riemann solver shares with the exact Riemann solution is that 


Zs p R p =F(U J )-F(U i ) 


so that (6.9) becomes 


( 6 . 10 ) 


~khij(F (U j) - F (Uj)). 


( 6 . 11 ) 
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This is the flux difference between C, and C j which has been split via (6.7) and (6.8). 

The method as described so far is thus equivalent to traditional flux-difference splitting methods. By 
integrating the conservation laws around each cell, it is easily shown that such methods are conservative for any 
splitting of the flux difference. A variety of such methods have been proposed, most of which share the 
property that flux differences are split only between C, and Cj with no contribution to other cells. This leads to 
a restriction on the time step, since clearly the factor multiplying R } in (6.7), for example, should be no greater 
than 1 for stability. If A, is very small, this gives a severe restriction on k . 

The advantage of the wave propagation interpretation of Figure 6.1 is twofold. First, it makes it clear that 
in some circumstances other cells should also be affected. More importantly, it makes it easy to determine 
exactly what that effect should be. Figure 6.2 shows a typical wave propagating over a different grid. Each cell 
value should clearly be updated by sgr\(s p )R p times the ratio of the shaded area within the cell to the total cell 
area. Note that this ratio is always less than or equal to 1. Splitting each wave in this manner and adding up all 
increments to all cells gives the same total increment (6.9) and hence the method remains conservative, while 
being much more stable in a situation such as Figure 6.2 than a method that increments only C,- and C ; . 

Instabilities can still arise, however, as illustrated by the following simple example. Consider the scalar 
linear advection equation 


u, + u x + Uy = 0 (6.12) 

on a uniform rectangular mesh with initial data consisting of a sawtooth discontinuity at 45° to the mesh: 


t/yH 


1 if i + j £ Af 
0 if t + y Af + 1 


for some Af . If k < hl2 then the method described above is stable and in fact for k - hi 2 gives the exact 
solution, a wave propagating through one additional diagonal of mesh cells in each time step: 


_ 1 if i + j < Af + 1 

^ jO ifi+y£Af+2 

For larger k, however, the method is unstable. For example, if k = h then the wave should propagate through 
two additional diagonals in each step. Instead the waves propagating in the x- and y -directions overlap and 
give the value Cfy = 2 for i + y = Af + 1 with no propagation into the second diagonal. In later time steps this 
instability grows exponentially. 


In order to avoid this difficulty, it appears to be necessary to include tangential wave propagation as well 
as normal propagation. For equation (6.12) it is easy to see how to include tangential propagation. All waves 
should propagate with tangential velocity 1 as well as normal velocity 1. Figure 6.3 shows this wave 
propagation in the case k = h 12, giving 

«■ 

1 if i + j < Af + 1 
0.75 ifi+y =Af + 1 
= ' 0.25 if i + y = Af + 2 
0 if i + j > Af + 2 
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The method now gives exact propagation with k = h and indeed, except for smearing introduced by the 
averaging process, this modified method gives exact propagation for any time step k . 

For nonlinear systems of equations the tangential propagation is introduced as follows. After solving the 
normal Riemann problem we have waves R p traveling with speeds s p normal to the cell face. Suppose now that 
we split up the wave R p into subwaves R p l \ R®, • • • , R p (m) with 

£R p (?) = R p (6.13) 

9=1 

and assign tangential velocities s p l) ‘ • • • , s p n \ The wave in Figure 6.2, for example, might be split into two 
waves /? p (1) and R p (2) with tangential speeds r p (1) and s p (2) as shown in Figure 6.4. Note that the area swept out 
by each wave is I s p \kh ijt the same as the area of the original wave in Figure 6.2, independent of the tangential 
speeds s^\ If each of these waves is propagated as before to increment neighboring mesh cells, the total 
increment will be 


-khijlSpR™ + s p R p ®] = -khijSpRp (6.14) 

in view of (6.13). Hence the method remains conservative for an arbitrary splitting /? p (<?) and arbitrary speeds 
s p (,) provided only that (6.13) is satisfied. 

A natural way to choose the splitting is by solving a tangential Riemann problem with initial data having 
a discontinuity of magnitude R p . For the equation we take the tangential portion of (6.3), 

u, + G («)„ = <). (6.15) 

Figure 6. 1 suggests a possible set of initial data. The wave R p propagates with the state 

Ui + 2 R q (6.16a) 

q<p 

to one side and 

Uj-^Rq (6.16b) 

q>p 

to the other. Taking these as the left and right states for the equation (6.15) gives a discontinuity of magnitude 

R p which will then be resolved into waves R p (1) , • • ■ , R p (m) propagating in the q direction with speeds 

„(D ... ,(«) 

^ • 

For the linear equation (6.12), m = 1, s j = 1 and s f 1 * = 1 for all i , j . This gives the correct propagation 
as described above. More generally, consider a linear system 


u, +Au x +BUy = 0. 

First suppose that A and B have identical eigenvectors and hence are simultaneously diagonalizable. Then 
each wave R p obtained by splitting Uj-Uj into eigenvectors of A is already an eigenvector of B as well, so 
that R p (p) -R p and R p (,) = 0 for q *p. The normal and tangential speeds s p and 5 p (p) are the corresponding 
eigenvalues of A and B respectively. In this case the method gives exact propagation (followed by smearing 
due to the averaging process) for any time step k and we are essentially using the method of characteristics. 
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If A and B arc not simultaneously diagonalizable, or for nonlinear systems, each wave R p will in general 
be split into m waves R p q) . Except in special circumstances the method will not give exact propagation and the 
time step will need to be restricted in order to obtain reasonable results. Also note that when waves from 
different interfaces overlap, they are combined by linear superposition as in the one dimensional algorithm. 
This incurs additional errors on nonlinear problems. A serious error analysis of the method has not yet been 
attempted, but the numerical results of Section 9 give some indication of its capabilities. 

In the special case of a uniform Cartesian grid, with =h and A, * A 2 , the wave propagation becomes 
relatively simple, particularly if the Courant number kh\s mn \lh 2 is less than 1. Then each wave affects at 
most two cells, as in Figure 6.3, and the areas of intersection are easy to compute. Although we don’t present 
the details, the method can be written as a nine-point finite difference method in conservation form. If 
tangential wave propagation is suppressed, the method reduces to a five-point scheme. 


7. High resolution in two dimensions. 

Greater accuracy for smooth flows and better resolution of discontinuities can be obtained by introducing 
piecewise linear approximations as in the one dimensional method. Ideally one might like to use bilinear 
approximations matching slopes in the x- and y -directions simultaneously, but it turns out that a simpler 
approach is sufficient to reproduce the Lax-Wendroff method for scalar problems on uniform grids. 
Introducing minmod slope limiting gives a stable second order correction that is easy to generalize and 
implement for nonlinear systems on arbitrary grids. 

First consider a uniform Cartesian grid with Courant number less than one and suppress tangential 
propagation. Then we obtain a five-point flux-splitting method as discussed above. We can think of 
implementing this by applying the first order one dimensional algorithm of Section 3 simultaneously in the x- 
and y -directions. It is natural then to also apply the second order correction waves introduced there. In 
handling the interface between C/ iy and U i+l j, for example, we obtain waves R u ■ • ■ ,R m propagating in the 
x -direction. Suppose that R p propagates with speed s p > 0. Then U M j is incremented by 



according to the first order algorithm. For the second order correction we propagate the wave u (2) from Section 
3, viewed now as a two dimensional wave that is constant in y over the grid cells in question and piecewise 

1 k 

linear in x. Thus U i+1 j is incremented by — — s p (h - ks p )G p while U it is incremented by 

^ fl 

1 k 

— — — s p (h -ks p )o p , where o p is the slope vector in the x -direction, again defined by one of (3.9). Wave 

4* fl 

propagation in the y -direction is handled in an analogous manner. 

If we choose the Lax-Wendioff slope (3.9b), o p =R p /h, then for scalar equations the method just 
described reduces to the following two dimensional version of the Lax-Wendroff method: 


u t = U„ - + g(U iJ+1 ) - giUij.O] 

k 2 \(nu M jynVij)? (/•(^>-m-i, / )) 2 

+ 2 h 2 [ Uij-Ui-Xj 
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(g(U lJ+ i)-g(U tJ )) 2 OK^HK^- 1)) 2 ' 

+ Uij+x-Uij + Uv-Utj . , / 

This method is second order accurate. As before, in order to eliminate oscillations it is better to use the 
minmod slope (3.9c) in practice, giving a method that remains second order accurate except at extrema. 

When tangential propagation is included, it is possible to split up the second order correction wave into 
subwaves and propagate these in the tangential direction as well. Here we adopt a simpler approach and always 
propagate the correction wave normal to the interface. On uniform grids the second order correction is then 
very easy to implement and adds little cost to the algorithm. The numerical results of Section 9 indicate that 
this is sufficient to give increased resolution of two-dimensional phenomena. 

On unstructured grids, a high resolution correction can be defined in various ways. Here is one untested 
possibility, motivated by the above treatment of uniform grids. Let § and t| be the local coordinate directions 
defined by the interface between Cj and C r For the pth wave, define a slope vector a p in the ^-direction 
(normal to the interface). The simplest choice would be the "Lax-Wendroff ' slope 

a p = Rp/h * (7.2) 

where h* is some measure of the normal distance between cells C, and Cj, such as the difference in %- 
coordinates of the centers of mass of these two cells. Then a piecewise linear wave with this slope in the in- 
direction, zero slope in the T|-direction, and total integral zero can be propagated in the ^-direction with velocity 
s p . The shape of the region where the slope is o p is still undetermined. Following the one dimensional 
method, we could choose the region to be the cell C, if s p > 0 and Cj if s p < 0. To be specific, assume that 
s p >0 and let (^.tlc) be the center of mass of the cell Cj, in the local coordinates determined by the interface. 
Then the high resolution correction, which must be averaged onto the cells, is the piecewise linear function 

... f£ - s p k - % c )c p if (5 - s p k ,r|) e C t 

u ” j0 otherwise. 

In problems with shocks one would again like to introduce a slope limiter instead of using (7.2) directly. 
The best way to do this on unstructured grids has not yet been determined. 

8. Boundary conditions in two dimensions. 

Nonreflecting outflow boundary conditions are implemented in two dimensions in exactly the same way 
as in one dimension. Waves propagating out of the computational domain at such a boundary simply disappear 
and no new waves are propagated inward from boundary interfaces. 

At a solid wall boundary, the proper boundary conditions for the Euler equations are obtained by 
specifying zero normal velocity. If the boundary is flat, this boundary condition can be implemented in the 
same way as was done in Section 4 for the one dimensional Euler equations. We will consider this simple case 
first and see how to reinterpret it in a manner more suited to arbitrary boundaries. For flat boundaries, the grid 
cells near the boundary can be reflected across the boundary to form fictitious cells as in Figure 8.1. We assign 
values of U in these exterior cells based on the values in the mirror-image interior cells so that the density, 
pressure, and tangential velocity in the exterior cell agree with those values in the interior while the normal 
velocity is negated. Then we can simply compute over the extended domain using the algorithm described 
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above, ignoring the boundary. Each wave from the interior which crosses the boundary will be matched by a 
wave from the exterior in such a way that the normal velocity in mirror-image cells adjacent to the boundary 
will remain negatives of one another, simulating zero normal velocity on the wall itself. 

Alternatively, rather than creating fictitious exterior cells and propagating waves from their interfaces to 
match interior waves, we can simply reflect the interior waves when they hit the boundary as suggested by 
Figure 8.1. If the wave carries increments R p then the reflected wave carries the same increments of density, 
pressure, and tangential velocity. The increment of normal velocity is negated. 

The latter approach has the advantage that it extends more easily to general boundaries. If the boundary 
is not flat, creating globally well-defined fictitious cells by reflecting the interior cells is not possible. It is, 
however, still possible to reflect waves at the boundary provided only that we assume the boundary is locally 
flat where it intersects the wave, as can be done to good approximation if the boundary is smooth. 

In addition to handling the reflection of waves arising from interior interfaces, it is also necessary to 
generate a wave from the cell face defining the boundary itself. We solve a one dimensional Riemann problem 
normal to the boundary. The left state comes from the interior cell. The right state is the same except that the 
normal velocity is negated, as motivated by Figure 8.1 and the previous discussion. Because the two states are 
related in this manner, the solution of the Riemann problem will give only one wave propagating into the 
computational domain, with zero normal velocity in its wake. Of the other three waves, two will have zero 
velocity and the third will propagate out of the computational domain and is ignored. 

9. Numerical results. 

To date the method has been implemented only on a Cartesian grid with "small" time steps, meaning the 
Courant number is less than one relative to the area h 1 of the regular grid cells. However, boundaries are 
allowed to cut through this grid and the area of the resulting boundary cells may be arbitrarily small. The 
boundary is approximated by piecewise linear segments connecting the points where the boundary crosses the 
grid lines. Consequently all cells are polygons with at most five sides, since we assume each cell contains at 
most one boundary segment. 

Because of the Courant number restriction, wave propagation in the interior cells is easily implemented 
as discussed in Section 6. It is only near the boundaries that waves may intersect several cells in an 
unpredictable manner. Each wave is also a polygon which can be represented by a list of its vertices. This 
polygon can be intersected with the polygonal grid cells and the area of intersection calculated using standard 
algorithms of computational geometry. A wave that strikes the boundary is reflected by intersecting it with a 
half space which locally approximates the exterior. This gives a polygon representing the portion of the wave 
exterior to the domain which is then reflected through the boundary of the half space to obtain an interior 
polygon. 

This basic algorithm has been fully implemented. The high resolution correction terms have been 
properly included only at the regular grid cells in the interior. At the boundary an approximation is used that 
gives reasonable results but can probably be improved upon. 

The test cases presented here involve only flat boundaries. The handling of boundary conditions is then 
equivalent to reflecting the interior cells across the boundary. This minimizes the possible negative effects of 
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the boundary conditions and gives a reasonable test of the basic algorithm in situations where there are irregular 
small grid cells. 

Allowing arbitrarily small boundary cells poses some numerical difficulties not associated with the 
theoretical stability properties of the method. If the size of a grid cell is on the order of the machine rounding 
error, then the area of the intersection of that cell with a wave cannot be accurately computed. This can lead to 
gross errors that will eventually cause noticeable errors in other grid cells and possibly nonphysical flow 
quantities. In order to avoid such problems, we have chosen to ignore any boundary grid cells with area less 
than I0~ 2 h 2 . This amounts to introducing imperceptible holes in the otherwise solid boundary. The 
experiments reported here were performed on a Vaxstation II in single precision. 

The first test case is a one dimensional shock tube problem in a channel at 40° to the grid. We take the 
same initial data as in Section 5. Figure 9.1 shows the boundary cells of the computational grid and contour 
lines of constant density for the numerical solution at t = 0.3, using the high resolution method with h - 1/40 
and timestep k = 0.0075. This gives a Courant number of roughly 0.5 relative to the regular cells, though it is 
nearly 100 times larger for the smallest irregular cells that are not ignored. Figure 9.1a shows the results 
obtained with the basic first order method. Figure 9.1b shows the results when the high resolution correction 
waves discussed in Section 7 are included. Both figures show contour lines of density as well as the irregular 
grid cells adjacent to the boundary. Figure 9.2 shows corresponding plots of the density along the lower wall of 
the channel. Values in each boundary cell are plotted vs. distance along the wall. No attempt has been made to 
extrapolate to the actual boundary location, or to account for the nonuniformity of these boundary cells. 
Consequently some nonsmoothness is visible in this boundary data (and even more so in Figure 9.4 below). 
Plots of grid values down the center of the channel would show much smoother results. 

Smearing is clearly reduced by introducing the high resolution correction waves. In both cases the 
contour plots show that the method is capable of maintaining one-dimensionality of the solution, with 
essentially the same accuracy in the irregular cells as in the interior. The small boundary cells cause no stability 
problems. 

The second test example is an unsteady shock reflection from a 27° ramp at Mach 2.03, with single Mach 
reflection. Experimental as well as numerical results for this problem can be found in the paper of Glaz, 
Collela, Glass, and Deschambault (Case 2 in [9]). The initial flow conditions are p = 0.387, v = 0 ,p =33.3 
ahead of the shock and p = 1.05, v = -14.06, p = 154.5 behind the shock, where v if the velocity normal to the 
shock. The shock was initially distance 0.05 from the ramp comer and results are presented at time 0.275 with 
h = 1/80 and k = 2.5X10 -4 (1 10 time steps). The results presented here were computed with the high resolution 
method. Computations with the first order version have also been done and give the same qualitative results 
with much more smearing. 

For comparison purposes, two different orientations of the ramp are considered. In Case I (Figure 9.3a), 
the grid is aligned with the ramp. This is the orientation often used in other codes to avoid irregular grid cells in 
the region of primary interest (e.g. [9], [29]). In Case II (Figure 9.3b) the grid is aligned with the incident shock 
and the ramp cuts through the grid obliquely. In each case 30 equally spaced density contours are plotted. The 
results agree closely and have essentially the same resolution. 

Figure 9.4 shows a plot of relative density (p/p 0 , where p 0 = 0.387) along the lower wall. The solid line 
shows the results from Case I while the + symbols show results from Case II. The squares are experimental 
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results taken from [9]. The results are in good agreement and also compare reasonably well with numerical 
results of [9]. The discontinuities obtained here are not as sharp, particularly the slip line connecting the triple 
point to the wall. In large part this is due to the coarseness of the grid (essentially 80x40 grid points compared 
with 350x130 in [9]). Figure 9.4 shows that even in Case II, where the ramp is oblique to the grid, there are at 
most two points in the shocks. 

10. Conclusions. 

An approach for deriving high resolution numerical methods has been presented that allows 
generalization to stable methods on arbitrary grids. Numerical results in two space dimensions show that the 
accuracy in irregular grid cells is comparable to what is obtained with regular cells. In particular, for Cartesian 
grids cut by irregular boundaries, it is possible to obtain stable high resolution results with a time step chosen 
based on the interior regular cells. 

Another application of these ideas is to mesh refinement using overlapping rotated Cartesian grids, as in 
the work of Berger[l,2], Again irregular cells are created at the grid interfaces and the present approach gives a 
way of deriving conservative interface conditions. This is currently being studied by the author and M. Berger. 

A similar approach toward achieving stability has been taken by Chem and Colella[5]. They give an 
algorithm for shock tracking across a Cartesian grid in which fluxes are restributed from small grid cells to 
neighboring grid cells to insure stability. Their redistribution procedure is rather ad hoc compared to the wave 
propagation approach taken here, but presumably has the advantage of being much cheaper. A major 
shortcoming of the present method is its expense. If tangential wave propagation is used in two dimensions, 
then at each cell interface m + 1 Riemann problems must be solved and m 2 waves propagated. Even using 
Roe’s approximate Riemann solver (which is easy to vectorize over multiple Riemann problems, particularly on 
a Cartesian grid), this may be prohibative for a practical algorithm. On the other hand, it does not seem that so 
much information should be necessary, and efforts are currently under way to find simplifications of the 
algorithm by using a smaller number of approximate waves that carry the essential information. One possibility 
is to use a different kind of wave decomposition in place of the Riemann solutions. Some of the ideas of 
Roe[24] or Deconinck, Hirsch and Peuteman[8] might be useful in this regard. 

Another practical possibility is to use the algorithm developed here only at irregular boundaries of a 
Cartesian grid, coupling it to a more efficient standard method in the regular interior cells. 

On the other hand, in some applications the approach discussed here may have advantages over standard 
finite volume or dimensionally split methods even on regular grids. The fact that waves are propagated oblique 
to the grid (when tangential propagation is included) may be advantageous in modeling complex two 
dimensional phenomena. Tests have shown that including tangential propagation, even when not required for 
stability, increases the accuracy with which waves oblique to the grid can be tracked. 

In some special applications it could also prove useful to use Courant numbers larger than one even on a 
regular grid. For example, in problems with widely varying wave speeds where the phenomena of primary 
interest occurs on the slow scale, the fast waves may carry little information but will limit the time step with 
traditional methods. Using the wave propagation approach allows one to choose a time step appropriate for the 
slow scale phenomena. This could lead to a simultaneous increase in effeiciency and reduction in the numerical 
dissipation of slow waves. 



- 18 - 


These possibilities are currently under study and will be reported on elsewhere. 
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Figure 2.1. The piecewise linear funciton at time t. 
at time f n+1 . Squares indicate the grid values {7y_ It • 
represents the contribution of this wave toward updating 


and after propagation, 
• . The shaded region 
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Figure 6 . 1 . Two grid cells C* and C ; - from an unstructured grid. Wave propaga- 
tion from the cell interface is determined by solving a one dimensional Riemann 
problem in the ^-direction, giving rise to waves propagating with velocities s lt s 2 > 
and S3. 



Figure 6.2. Propagation of a wave over a general unstructured grid. Cell 
values are updated in every cell the wave intersects. The contribution to each cell 
is proportional to the ratio of the shaded area within that cell to the total cell area. 
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boundary 



. Figure 8.1. Cartesian mesh cut by a flat boundary. Reflecting the interior cells 
through this boundary gives the fictitious exterior ceils. Each wave propagating 
across the boundary from an interior cell interface is matched by a wave propagat- 
ing into the region from a fictitious interface. This can also be viewed as a 
reflection of the interior wave at the boundary. 










Figure 9.3 Density contours for Mach 2.03 shock reflection off a 27° ramp, (a) 
Case I: grid aligned with ramp, (b) Case II: grid aligned with incident shock. 



Figure 9.4 Density along the lower boundary. The solid line shows results from 
Case I, the + symbols are results from Case II. The squares are experimental 
results taken from [9]. 








